#######################
#######################
###Main Paper Models###
#######################
#######################
library(readstata13)
library(countrycode)
library(foreign)
library(dplyr)
library(doBy)
library(DataCombine)
library(AER)
library(lfe)
library(lmtest)
library(survival)
library(sandwich)
library(stargazer)
library(ggplot2)
library(data.table)
#Set random number generator values
set.seed(123123)
##Read data
#Set directory
setwd("~/OneDrive - Indiana University/FromGoogle/DiseaseFood/AJARESubmissionCode_4_24/")
#Read dataset
dis.dat <- read.csv("diseasefooddatfull12822.csv")

#############################
###Create a time indicator###
#############################
#Create a time indiactor for lengh
t <- 0:(length(unique(dis.dat$year)))
#Match with years
year <- 1995:2020
#Match by year
t2 <- data.frame(cbind(year,t))
summary(t2)
#Redefine name
names(t2) <- c("year", "t")
#Join the new indicator into the data
dis.dat <- left_join(dis.dat, t2)
summary(dis.dat)

###################################################
###Create logged versions of relevant indicators###
###################################################
dis.dat$logpopulation <- log(dis.dat$population+1)
dis.dat$logGDPpc <- log(dis.dat$GDPpc+1)

#####################################################
###Create 1000000s versions of relevant indicators###
#####################################################
dis.dat$comb.fert.use_tonnes1000000 <- (dis.dat$comb.fert.use_tonnes)/1000000
dis.dat$comb.fert.imp_tonnes1000000 <- (dis.dat$comb.fert.imp_tonnes)/1000000
dis.dat$nitrogen.use_tonnes1000000 <- (dis.dat$nitrogen.use_tonnes)/1000000
dis.dat$nitrogen.imp_tonnes1000000 <- (dis.dat$nitrogen.imp_tonnes)/1000000
dis.dat$phosphate.use_tonnes1000000 <- (dis.dat$phosphate.use_tonnes)/1000000
dis.dat$phosphate.imp_tonnes1000000 <- (dis.dat$phosphate.imp_tonnes)/1000000
dis.dat$cattle_meat_tonnes1000000 <- (dis.dat$cattle_meat_tonnes)/1000000
dis.dat$chicken_meat_tonnes1000000 <- (dis.dat$chicken_meat_tonnes)/1000000
dis.dat$pigs_meat_tonnes1000000 <- (dis.dat$pigs_meat_tonnes)/1000000
#Create a combined fertilizer indicator for phosphates and nitrates
dis.dat$comb.np.use_tonnes <- dis.dat$nitrogen.use_tonnes+dis.dat$phosphate.use_tonnes
dis.dat$comb.np.imp_tonnes <- dis.dat$nitrogen.imp_tonnes+dis.dat$phosphate.imp_tonnes
#Create a 1000000s version
dis.dat$comb.np.use_tonnes1000000 <- dis.dat$nitrogen.use_tonnes1000000+dis.dat$phosphate.use_tonnes1000000 
dis.dat$comb.np.imp_tonnes1000000 <- dis.dat$nitrogen.imp_tonnes1000000+dis.dat$phosphate.imp_tonnes1000000
##Create pathogen group indicators
#Flus
dis.dat$flus_all <- dis.dat$H1N1+dis.dat$H5N1+dis.dat$H5N6+dis.dat$H7N9+dis.dat$flu_other 
dis.dat$lagflus_all <- dis.dat$lagH1N1+dis.dat$lagH5N1+dis.dat$lagH5N6+dis.dat$lagH7N9+dis.dat$lagflu_other 
#Fevers
dis.dat$fevers_all <- dis.dat$chikungunya+dis.dat$dengue+dis.dat$ebola_reston+dis.dat$y_fever+dis.dat$enterovirus+dis.dat$j_enchaphalitis+dis.dat$nipah+dis.dat$zika 
dis.dat$lagfevers_all <- dis.dat$lagchikungunya+dis.dat$lagdengue+dis.dat$lagebola_reston+dis.dat$lagy_fever+dis.dat$lagenterovirus+dis.dat$lagj_enchaphalitis+dis.dat$lagnipah+dis.dat$lagzika
#Respiratory syndrome
dis.dat$resp_all <- dis.dat$MERS+dis.dat$SARS
dis.dat$lagresp_all <- dis.dat$lagMERS+dis.dat$lagSARS

##Subset relevant years only
dis.dat <- subset(dis.dat, year>1995&year<2020)

################
################
###OLS Models###
################
################

###############
####Baseline###
############### 
##Cattle
lm1.cat <- lm(outbreak_event~cattle_meat_tonnes1000000+
                t+factor(ccode), data=dis.dat)
#Cluster SEs
lm1.cat.sum <- coeftest(lm1.cat, vcov = vcovCL(lm1.cat, cluster = ~ccode))

##Chicken
lm1.chick <- lm(outbreak_event~chicken_meat_tonnes1000000+
                  t+factor(ccode), data=dis.dat)
#Cluster SEs
lm1.chick.sum <- coeftest(lm1.chick, vcov = vcovCL(lm1.chick, cluster = ~ccode))

##Pigs
lm1.pigs <- lm(outbreak_event~pigs_meat_tonnes1000000+
                 t+factor(ccode), data=dis.dat)
#Cluster SEs
lm1.pigs.sum <- coeftest(lm1.pigs, vcov = vcovCL(lm1.pigs, cluster = ~ccode))


###########
####Full###
###########
##Cattle
lm2.cat <- lm(outbreak_event~cattle_meat_tonnes1000000+lagoutbreak_event+
                logpopulation+logGDPpc+
                t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.cat.sum <- coeftest(lm2.cat, vcov = vcovCL(lm2.cat, cluster = ~ccode))

##Chicken
lm2.chick <- lm(outbreak_event~chicken_meat_tonnes1000000+lagoutbreak_event+
                  logpopulation+logGDPpc+
                  t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.chick.sum <- coeftest(lm2.chick, vcov = vcovCL(lm2.chick, cluster = ~ccode))

##Pigs
lm2.pigs <- lm(outbreak_event~pigs_meat_tonnes1000000+lagoutbreak_event+
                 logpopulation+logGDPpc+
                 t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.pigs.sum <- coeftest(lm2.pigs, vcov = vcovCL(lm2.pigs, cluster = ~ccode))

##Export OLS with CSEs table to Latex
stargazer(lm1.cat.sum, lm1.chick.sum, lm1.pigs.sum, lm2.cat.sum, lm2.chick.sum, lm2.pigs.sum)
#Obs and diagnostics
stargazer(lm1.cat, lm1.chick, lm1.pigs, lm2.cat, lm2.chick, lm2.pigs)

###############
###############
###IV Models###
###############
###############

###############
####Baseline###
###############
####Cattle
iv1.cat <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+t+factor(ccode)|
                      comb.np.use_tonnes1000000+t+factor(ccode), data=dis.dat)
summary(iv1.cat, diagnostics=T)
#Cluster results
iv1.cat.sum <- coeftest(iv1.cat, vcov = vcovCL(iv1.cat, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv1.cat.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv1.chick <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+t+factor(ccode)|
                    comb.np.use_tonnes1000000+t+factor(ccode), data=dis.dat)
summary(iv1.chick, diagnostics=T)
#Cluster results
iv1.chick.sum <- coeftest(iv1.chick, vcov = vcovCL(iv1.chick, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv1.chick.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv1.pigs <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+t+factor(ccode)|
                      comb.np.use_tonnes1000000+t+factor(ccode), data=dis.dat)
summary(iv1.pigs, diagnostics=T)
#Cluster results
iv1.pigs.sum <- coeftest(iv1.pigs, vcov = vcovCL(iv1.pigs, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv1.pigs.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


###########
####Full###
###########
####Cattle
iv2.cat <-  ivreg(outbreak_event~cattle_meat_tonnes1000000+t+
                    lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                    comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat, diagnostics=T)
#Cluster results
iv2.cat.sum <- coeftest(iv2.cat, vcov = vcovCL(iv2.cat, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick <-  ivreg(outbreak_event~chicken_meat_tonnes1000000+t+
                      lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                      comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick, diagnostics=T)
#Cluster results
iv2.chick.sum <- coeftest(iv2.chick, vcov = vcovCL(iv2.chick, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs <-  ivreg(outbreak_event~pigs_meat_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode)|
                     comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs, diagnostics=T)
#Cluster results
iv2.pigs.sum <- coeftest(iv2.pigs, vcov = vcovCL(iv2.pigs, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("outbreak_event","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv1.cat.sum, iv1.chick.sum, iv1.pigs.sum, iv2.cat.sum, iv2.chick.sum, iv2.pigs.sum)
#Obs and diagnostics
stargazer(iv1.cat, iv1.chick, iv1.pigs, iv2.cat, iv2.chick, iv2.pigs)
#Corrected weak instruments
iv.main <- data.frame(cbind(iv1.cat.sum.stock,iv1.chick.sum.stock,iv1.pigs.sum.stock,iv2.cat.sum.stock,iv2.chick.sum.stock,iv2.pigs.sum.stock))
colnames(iv.main) <- c("Cattle (Base)", "Chicken (Base)", "Pigs (Base)", "Cattle (Full)", "Chicken (Full)", "Pigs (Full)")
iv.main



#########################################
#########################################
###Disaggregate by Disease Type -- OLS###
#########################################
#########################################

###########
####Flus###
###########
##Cattle
lm2.cat.flus <- lm(flus_all~cattle_meat_tonnes1000000+lagflus_all+
                     logpopulation+logGDPpc+
                     t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.cat.flus.sum <- coeftest(lm2.cat.flus, vcov = vcovCL(lm2.cat.flus, cluster = ~ccode))

##Chicken
lm2.chick.flus <- lm(flus_all~chicken_meat_tonnes1000000+lagflus_all+
                       logpopulation+logGDPpc+
                       t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.chick.flus.sum <- coeftest(lm2.chick.flus, vcov = vcovCL(lm2.chick.flus, cluster = ~ccode))

##Pigs
lm2.pigs.flus <- lm(flus_all~pigs_meat_tonnes1000000+lagflus_all+
                      logpopulation+logGDPpc+
                      t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.pigs.flus.sum <- coeftest(lm2.pigs.flus, vcov = vcovCL(lm2.pigs.flus, cluster = ~ccode))


############
###Fevers###
############
##Cattle
lm2.cat.fevers <- lm(fevers_all~cattle_meat_tonnes1000000+lagfevers_all+
                       logpopulation+logGDPpc+
                       t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.cat.fevers.sum <- coeftest(lm2.cat.fevers, vcov = vcovCL(lm2.cat.fevers, cluster = ~ccode))

##Chicken
lm2.chick.fevers <- lm(fevers_all~chicken_meat_tonnes1000000+lagfevers_all+
                         logpopulation+logGDPpc+
                         t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.chick.fevers.sum <- coeftest(lm2.chick.fevers, vcov = vcovCL(lm2.chick.fevers, cluster = ~ccode))

##Pigs
lm2.pigs.fevers <- lm(fevers_all~pigs_meat_tonnes1000000+lagfevers_all+
                        logpopulation+logGDPpc+
                        t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.pigs.fevers.sum <- coeftest(lm2.pigs.fevers, vcov = vcovCL(lm2.pigs.fevers, cluster = ~ccode))

##########################
###Respiratory syndrome###
##########################
##Cattle
lm2.cat.resp <- lm(resp_all~cattle_meat_tonnes1000000+lagresp_all+
                     logpopulation+logGDPpc+
                     t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.cat.resp.sum <- coeftest(lm2.cat.resp, vcov = vcovCL(lm2.cat.resp, cluster = ~ccode))

##Chicken
lm2.chick.resp <- lm(resp_all~chicken_meat_tonnes1000000+lagresp_all+
                       logpopulation+logGDPpc+
                       t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.chick.resp.sum <- coeftest(lm2.chick.resp, vcov = vcovCL(lm2.chick.resp, cluster = ~ccode))

##Pigs
lm2.pigs.resp <- lm(resp_all~pigs_meat_tonnes1000000+lagresp_all+
                      logpopulation+logGDPpc+
                      t+factor(ccode), data=dis.dat)
#Cluster SEs
lm2.pigs.resp.sum <- coeftest(lm2.pigs.resp, vcov = vcovCL(lm2.pigs.resp, cluster = ~ccode))

##Export OLS CSEs table to Latex
stargazer(lm2.cat.flus.sum, lm2.chick.flus.sum, lm2.pigs.flus.sum, lm2.cat.fevers.sum, lm2.chick.fevers.sum, lm2.pigs.fevers.sum, lm2.cat.resp.sum, lm2.chick.resp.sum, lm2.pigs.resp.sum)
#Obs and diagnostics
stargazer(lm2.cat.flus, lm2.chick.flus, lm2.pigs.flus, lm2.cat.fevers, lm2.chick.fevers, lm2.pigs.fevers, lm2.cat.resp, lm2.chick.resp, lm2.pigs.resp)


#########################################
#########################################
###Disaggregate by Disease Type -- IVs###
#########################################
#########################################

##########
###Flus###
##########
####Cattle
iv2.cat.flus <-  ivreg(flus_all~cattle_meat_tonnes1000000+t+
                    lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                    comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.flus, diagnostics=T)
#Cluster results
iv2.cat.flus.sum <- coeftest(iv2.cat.flus, vcov = vcovCL(iv2.cat.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.flus.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.flus <-  ivreg(flus_all~chicken_meat_tonnes1000000+t+
                      lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                      comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.flus, diagnostics=T)
#Cluster results
iv2.chick.flus.sum <- coeftest(iv2.chick.flus, vcov = vcovCL(iv2.chick.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.flus.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.flus <-  ivreg(flus_all~pigs_meat_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode)|
                     comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.flus, diagnostics=T)
#Cluster results
iv2.pigs.flus.sum <- coeftest(iv2.pigs.flus, vcov = vcovCL(iv2.pigs.flus, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("flus_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagflus_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagflus_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.flus.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

############
###Fevers###
############
####Cattle
iv2.cat.fevers <-  ivreg(fevers_all~cattle_meat_tonnes1000000+t+
                         lagfevers_all+logpopulation+logGDPpc+factor(ccode)|
                         comb.np.use_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.fevers, diagnostics=T)
#Cluster results
iv2.cat.fevers.sum <- coeftest(iv2.cat.fevers, vcov = vcovCL(iv2.cat.fevers, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("fevers_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagfevers_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.fevers.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.fevers <-  ivreg(fevers_all~chicken_meat_tonnes1000000+t+
                           lagfevers_all+logpopulation+logGDPpc+factor(ccode)|
                           comb.np.use_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.fevers, diagnostics=T)
#Cluster results
iv2.chick.fevers.sum <- coeftest(iv2.chick.fevers, vcov = vcovCL(iv2.chick.fevers, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("fevers_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagfevers_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.fevers.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.fevers <-  ivreg(fevers_all~pigs_meat_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode)|
                          comb.np.use_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.fevers, diagnostics=T)
#Cluster results
iv2.pigs.fevers.sum <- coeftest(iv2.pigs.fevers, vcov = vcovCL(iv2.pigs.fevers, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("fevers_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagfevers_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagfevers_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.fevers.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

#########################
###Respiratory syndrom###
#########################
####Cattle
iv2.cat.resp <-  ivreg(resp_all~cattle_meat_tonnes1000000+t+
                           lagresp_all+logpopulation+logGDPpc+factor(ccode)|
                           comb.np.use_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.cat.resp, diagnostics=T)
#Cluster results
iv2.cat.resp.sum <- coeftest(iv2.cat.resp, vcov = vcovCL(iv2.cat.resp, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("resp_all","cattle_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagresp_all","logpopulation","logGDPpc")])
fs = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(cattle_meat_tonnes1000000 ~ t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.cat.resp.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Chicken
iv2.chick.resp <-  ivreg(resp_all~chicken_meat_tonnes1000000+t+
                             lagresp_all+logpopulation+logGDPpc+factor(ccode)|
                             comb.np.use_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.chick.resp, diagnostics=T)
#Cluster results
iv2.chick.resp.sum <- coeftest(iv2.chick.resp, vcov = vcovCL(iv2.chick.resp, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("resp_all","chicken_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagresp_all","logpopulation","logGDPpc")])
fs = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(chicken_meat_tonnes1000000 ~ t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.chick.resp.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]

####Pigs
iv2.pigs.resp <-  ivreg(resp_all~pigs_meat_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode)|
                            comb.np.use_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data=dis.dat)
summary(iv2.pigs.resp, diagnostics=T)
#Cluster results
iv2.pigs.resp.sum <- coeftest(iv2.pigs.resp, vcov = vcovCL(iv2.pigs.resp, cluster = ~ccode))
##Testing  weak instruments, adjusting for heteroskedasticity
ttt <- na.omit(dis.dat[,c("resp_all","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagresp_all","logpopulation","logGDPpc")])
fs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
# null first-stage (i.e. exclude IVs):
fn = lm(pigs_meat_tonnes1000000 ~ t+lagresp_all+logpopulation+logGDPpc+factor(ccode), data = ttt)
lmtest::waldtest(fs, fn)$F[2]
# F-test robust to heteroskedasticity
iv2.pigs.resp.sum.stock <- lmtest::waldtest(fs, fn, vcov = vcovHC(fs, type="HC0"))$F[2]


##Export IV models table to Latex
#Estimates
stargazer(iv2.cat.flus.sum, iv2.chick.flus.sum, iv2.pigs.flus.sum, iv2.cat.fevers.sum, iv2.chick.fevers.sum, iv2.pigs.fevers.sum, iv2.cat.resp.sum, iv2.chick.resp.sum, iv2.pigs.resp.sum)
#Obs and diagnostics
stargazer(iv2.cat.flus, iv2.chick.flus, iv2.pigs.flus, iv2.cat.fevers, iv2.chick.fevers, iv2.pigs.fevers, iv2.cat.resp, iv2.chick.resp, iv2.pigs.resp)
#Corrected weak instruments
iv.main.flus <- data.frame(cbind(iv2.cat.flus.sum.stock,iv2.chick.flus.sum.stock,iv2.pigs.flus.sum.stock, iv2.cat.fevers.sum.stock,iv2.chick.fevers.sum.stock,iv2.pigs.fevers.sum.stock, iv2.cat.resp.sum.stock,iv2.chick.resp.sum.stock,iv2.pigs.resp.sum.stock))
colnames(iv.main.flus) <- c("Cattle (Flu)", "Chicken (Flu)", "Pigs (Flu)","Cattle (Fever)", "Chicken (Fever)", "Pigs (Fever)", "Cattle (Resp.)", "Chicken (Resp.)", "Pigs (Resp.)")
iv.main.flus




                      
                      